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Abstract 

This paper focuses on the oscillation threshold of single reed 
instruments. Several characteristics such as blowing pres- 
sure at threshold, regime selection, and playing frequency 
are known to change radically when taking into account the 
reed dynamics and the flow induced by the reed motion. 
Previous works have shown interesting tendencies, using an- 
alytical expressions with simplified models. In the present 
study a more elaborated physical model is considered. The 
influence of several parameters, depending on the reed prop- 
erties, the design of the instrument or the control operated 
by the player, are studied. Previous results on the influence 
of the reed resonance frequency are confirmed. New results 
concerning the simultaneous influence of two model param- 
eters on oscillation threshold, regime selection and playing 
frequency are presented and discussed. The authors use a 
numerical continuation approach. Numerical continuation 
consists in following a given solution of a set of equations 
when a parameter varies. Considering the instrument as a 
dynamical system, the oscillation threshold problem is for- 
mulated as a path following of Hopf bifurcations, generalis- 
ing the usual approach of the characteristic equation, as used 
in previous works. The proposed numerical approach proves 
to be useful for the study of musical instruments. It is com- 
plementary to analytical analysis and direct time-domain or 
frequency-domain simulations since it allows to derive infor- 
mation that is hardly reachable through simulation, without 
the approximations needed for analytical approach. 
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1 Introduction 

Woodwind instruments have been intensively studied since 
Helmholtz in the late 19*^ century. In particular, the 
question of the oscillation threshold has been the focus of 
many theoretical, analytical, numerical and experimental 
studies[l-7j. Following the works of Backus and Benade, 
Worman[8j founded the basis of the small oscillation theory 
for clarinet-like systems, and Wilson and Beavers [9j (W&B) 
greatly extended his work on the oscillation threshold. 

Thompson[10j mentioned the existence of an additional 
flow, related to the reed motion, and the importance of the 
reed resonance on intonation as early as 1979, but only re- 
cent literature such as Kergomard and Chaigne[llj (chap. 
9), Silva et al.[12j and Silva[13j extended the previous works 
of Wilson and Beavers to more complex models, taking into 
account the reed dynamics and the reed-induced flow. 

Whereas the static regime of the clarinet (i.e. when 
the player blows the instrument but no note is played) 
can be derived analytically the stability of this regime 
and the blowing pressure threshold at which a stable pe- 
riodic solution occurs is not always possible to write an- 
alytically depending on the complexity of the model un- 
der consideration. In the latter two references, authors 
analyse the linear stability of the static regime through a 
characteristic equation. It is a perturbation method based 
on the equations written in the frequency domain, that 
suppose an arbitrarily small oscillating perturbation (with 
unknown angular frequency u) around the static solution. 
Substituting the linear expressions of the passive compo- 
nents - involving the acoustical impedance Z defined as: 
P{uj)=Z{uj)U{uj) and the reed's mechanical transfer func- 
tion D defined as: X{(j)=D{uj)P{uj), where P{uj),U{uj) are 
the acoustical pressure and volume flow at the input end of 
the resonator and X{uj) is the reed tip displacement - into 
the coupling equation linearised around the static solution 
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{U = /^*^|^^^^^^ (X, P)) and balancing the first harmonic 
of the oscillating perturbation, one gets the characteristic 
equation that the angular frequency uj and the blowing pres- 
sure 7 must satisfy at the instability threshold. 

As pointed out by Silva[13j, the convergence towards an 
oscillating regime with angular frequency uj after destabil- 
isation of the static regime is not certain and requires to 
consider the whole non linear system and compute periodic 
solutions beyond the threshold, in order to determine the 
nature of the Hopf bifurcation (direct or inverse). However, 
in the case of a direct bifurcation, the oscillation threshold 
(for a slowly increasing blowing pressure) corresponds to the 
solution of the characteristic equation with the lowest blow- 
ing pressure jth- Periodic oscillations emerge at this point 
with an angular frequency ujth that depends only on the 
other model parameters values. 

In the present paper, we propose a different method for 
computing the oscillation threshold that proves to be faster, 
more general and more robust: numerical continuation. 
Numerical continuation consists in following a given solu- 
tion of a set of equations when a parameter varies. To 
the authors knowledge, it is the first time this numerical 
method is applied to the physics of musical instruments. 
While time-domain or frequency-domain simulations allow 
to study complex physical models with at most one vary- 
ing parameter, the method proposed in this work allows to 
investigate the simultaneous influence of several parameters 
on the oscillation threshold. 

The general principle of the method is the following : we 
first follow the static regime while increasing the blowing 
pressure, all other parameters being kept constant, and de- 
tect the Hopf bifurcations. Then, we follow each Hopf bi- 
furcation when a second parameter fi is allowed to vary (/i 
being, for instance, the reed opening parameter or the reed 
resonance frequency). Finally, the resulting branch of each 
bifurcation is plotted in the plane {ja^jth) for determination 
of the oscillation threshold and of the selected regime, and 
in the plane (/i, Uth) to determine the corresponding playing 
frequency. 

The paper is structured as follows. In section [2| a physical 
model of a clarinet is reviewed. The method is described in 
section [3) defining the continuation of static solutions and 
Hopf bifurcations. In section |4j the results are reviewed. 
First, basic results on the reed-bore interaction are com- 
pared with the previous works (W&B[9j and Silva et al.[12j). 
Then, extended results are presented, showing how the con- 
trol of the player and the design of the maker influence the 
ease of play and the intonation. 

2 Physical model of single reed in- 
struments 

The model used in this study is similar to the one used in 
Silva et al.[12j. It is a three-equation physical model that 
embeds the dynamical behaviour of the single reed, the lin- 



ear acoustics of the resonator, and the nonlinear coupling 
due to the air jet in the reed channel. It assumes that the 
mouth of the player (together with the vocal tract and lungs) 
provides an ideal pressure source, thus ignoring acoustic be- 
haviour upstream from the reed. 

2.1 Dynamics of the reed 

The reed is indeed a three-dimensional object. However, due 
to its very shallow shape, it can be reduced to a 2D-plate 
or even, considering the constant width, a ID-beam with 
varying thickness. Detailed studies of such models have been 
conducted (see for instance Avanzini and van Walstiin|14] 
for a distributed ID-beam model) and lumped models have 
been proposed. Mainly, two lumped models are widely used: 

• a one degree of freedom mass-spring-damper system, 
accounting for the first bending mode of the reed, 

• a simpler one degree of freedom spring, assuming a 
quasi-static behaviour of the reed. 

Most papers studying clarinet-like systems use the second 
model. This kind of simplification is made assuming that 
the reed first modal frequency is far beyond the playing fre- 
quency, thus acting as a low-pass filter excited in the low 
frequency range, almost in phase. However, some authors 
(see for instance W&B[9j or Silva[13j) showed that the reed 
dynamics, even rendered with a simple mass-spring-damper 
oscillator, cannot be ignored because of its incidence on the 
physical behaviour of the dynamical system. 

Thus, in the present study, the reed dynamics is rendered 
through a lumped, one-mode model. It reflects the first 
bending mode of a beam-like model with a reed modal an- 
gular frequency cj^, a modal damping g^, and an effective 
stiffness per unit area Ka. 

The reed is driven by the alternating pressure difference 
AP=Pm — P(t) across the reed, where Pm is the pressure 
inside the mouth of the player and P{t) is the pressure in the 
mouthpiece. Given the rest position of the tip of the reed 
Hq and the limit h = when the reed channel is closed, the 
tip of the reed position h{t) satisfies the following equation: 

^h'\t) = ho- h(t) - ^h\t) - ^ (1) 

Contact with the mouthpiece could also be modelled, us- 
ing a regularised contact force function (as presented by the 
authors in a conference paper[15j). It will not be discussed 
or used in this paper, as large amplitude oscillations will 
not be discussed. However, it is important to note that such 
phenomenon cannot be ignored when the whole dynamic 
range of the system is under study. 

2.2 Acoustics of the resonator 

The acoustic behaviour of the resonator of a clarinet is here 
assumed to be linear. It is modelled through the input 
impedance Zin=P/U^ where P is the acoustical pressure 
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and U the acoustical volume flow at the reed end of the res- 
onator. To keep a general formulation, we use a (truncated) 
complex modal decomposition of the input impedance (as 
described by Silva[13j) of the form: 



^in(^) 



c* 



(2) 



This leads, in the time domain, to a system of differen- 
tial equations, each governing a complex mode of pressure 
at the input of the resonator : 



P;^{t)=CnZ,U{t)^SnPn{t), 



(3) 



where Zc = pc/S is the characteristic impedance of the 
cylindrical resonator, whose cross-sectional area is S. The 
total acoustic pressure is then given by : 



Nrr 



P(i) = 2^Re(P„(t)). 



(4) 



In practice, the modal coefficients (Cn^Sn) can be com- 
puted in order to fit any analytical or measured impedance 
spectrum to the desired degree of accuracy (depending on 
the total number of modes Nm)- 

In the present study, we use an analytical formulation for a 
cylinder of length L = 57cm (if not stated otherwise), radius 
r = 7mm, taking into account the radiation at the output 
end and the thermoviscous losses due to wall friction, as 
given by Silva[13j. The effect of tone holes is not considered 
in this work. 

2.3 Nonlinear coupling 

As described by Hirschberg|16| and further confirmed exper- 
imentally by Dalmont et al.[T7] and Almeida et al.[T8], the 
flow inside the reed channel forms a jet that is dissipated by 
turbulence in the larger part of the mouthpiece leading to 
the following nonlinear equation: 



Uj{t) = sign{AP)Wh{t)i 



l2\AP\ 



(5) 



where W is the width of the reed channel, assumed to be 
constant, and Uj has the sign of AP. According to this 
model, the direction of the flow depends on the sign of the 
pressure drop across the reed AP. 

2.4 Reed motion induced flow 

Moreover, when an oscillation occurs (a note being played), 
the reed periodic motion induces an acoustic volume flow 
in addition to the one of the jet, as early described by 
Thompson|lQ) and lately studied by Silva et al.[T2]. Writing 
Sr the reed effective area Sr, the resulting flow Ur reads: 



(6) 



Notice the negative sign, due to the chosen representation 
where the reed tip position h is positive at rest and null 
when closing the reed channel. 

Some authors [191 [11] use another notation, the equivalent 
length correction AL, related to Sr as: 



AL 



We prefer to stick with the notation Sr but we will give the 
corresponding values of AL for comparison, when necessary. 

2.5 Global flow 

The global flow entering the resonator then reads: 



U{t) = sign{AP)Wh{t)A 



/2|AP| 



dh 



(7) 



2.6 Dimensionless model 



The above- described multi-physics model involves several 
variables. Expressed in SI units, their respective values are 
of very different orders of magnitude. However, as we will 
use numerical tools to solve the algebro-differential system 
composed of eq. [T] |3j |4] and [7| it is a useful precaution to 
use a dimensionless model. 

In the dimensionless model, each variable should be di- 
vided by a typical value, in order to scale the range of pos- 
sible values for that variable as close as possible to [0, 1]. 

Choosing the static difference of pressure necessary to 
close the reed channel Pm = Kah^ to scale all pressures, 
and the reed tip opening at rest to scale the reed tip 
position, the following new dimensionless variables and pa- 
rameter are defined: 

x{t) = h{t)/ho 
y{t) = h\t)/vo 
Pn{t)=Pn{t)/PM for n = l,2,...7V^ 

p{t) = P(t)/PM 
u{t) = U(t)ZjPM 

where '^o = h^ojr is a characteristic velocity of the reed. 

The system of algebraic and differential equations 
( 1|3|4|7 ) can be rewritten as follows (time-dependence is 
omitted, for a more concise writing): 



-y' 

Pn 
P 
U 



y 

l-X^p-J-QrV 

CnU^SnPn for n = 1 , 2, . . . A/"^ 



Nrr 



(8) 



2 E MVn) 



Slg] 



n(7-p)C^\/|7-p| 



Pm 



SrV, 



where 7 = Pm/PM is the dimensionless blowing pressure 
and C = ZcW ^2hQ/ pKa is the dimensionless reed open- 
ing parameter. The reader will notice that the parameter 
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related to the maximum flow through the reed channel, 
mainly depends on the geometry of the mouthpiece, the reed 
mechanical properties, as well as the player's lip force and 
position on the reed that control the opening. 

In the following sections we will compute the minimal 
mouth pressure ^th necessary to initiate self-sustained oscil- 
lations from the static regime and show how this threshold 
is modified when a second parameter of the model varies. 
When players change the way they put the mouthpiece in 
their mouth, several parameters of our model are assumed to 
vary : the control parameter the reed modal damping g^, 
its modal angular frequency cj^, as well as the effective area 
Sr of the reed participating to the additional flow. We will 
first investigate the influence of the dimensionless parame- 
ter krL = uOj-L/c^ comparing with previous works. Then the 
influence of <7r and Sr will also be investigated. 

3 Methods : Theoretical principles 
and numerical tools 

In this section, the method used is briefly reviewed from 
a theoretical point of view. Practical application of the 
method will be carried out in section |4| 

3.1 Branch of static solutions 

We consider an autonomous nonlinear dynamical system of 
the form: 

u'(t)=F(u(t),A) (9) 

where \i{t) is the state vector and A is a chosen parameter 
of the system equations. 

Let uq G be a fixed point of this system for Aq: 
F(uo, Ao) = 0. Under the hypothesis that uq is a regular so- 
lution, there exists a unique function u(A) that is solution of 
the previous equation for all A close to Aq, with u(Ao) = uq. 
(See Doedel[20j for formal definitions, as well as for more de- 
tails about regular solutions and what is meant by "close".) 

In other words, for A close to Aq, there is a continuum 
of equilibrium points that passes by uq and the branch of 
solutions is represented by the graph Ui = /(A), where Ui is 
a component of the state vector u. 

3.2 Continuation of static solutions 

Considering the differential equations in system (|8|, it can 
be rewritten in the form = F(u, A) where u is the state 
vector and A is one of the equations parameters (for instance 
the blowing pressure 7). Thus, a static solution is given 
by the algebraic system: F(u, A) = 0, where the algebraic 
equations of system (|8| can now be included. 

The branch of static solutions is computed numeri- 
cally using a classical path following technique based on 
Keller's pseudo-arc length continuation algorithm[21J as im- 
plemented in the softwares AUTO (Doedel and Oldeman 



[22J) and MANLAB (Karkar et al. [23j), both freely avail- 
able online. For comprehensive details about continuation 
of static solutions using these two numerical tools, please 
refer to the conference article [15j. 

3.3 Hopf Bifurcation 

Let us consider a static solution u of system (|9|. The sta- 
bility of this equilibrium point is given by eigenvalues of the 
jacobian matrix of the system: 



. If all the eigenvalues of F^ have a strictly negative real 
part, then the equilibrium is stable. If any of its eigenval- 
ues has a strictly positive real part, then the equilibrium is 
unstable. 

Several scenarios of loss of stability, along the branch u = 
/(A) are possible. One of them is the following : a unique 
pair of complex conjugate eigenvalues crosses the imaginary 
axis at {iuJth, -ioJth) for A = A^^. 

In that case, the system is said to undergo a Hopf Bifur- 
cation. It means that a family of periodic orbits starts from 
this point u(At^), with an angular frequency ujth- As for 
the clarinet model, choosing the blowing pressure 7 as the 
varying parameter, the first Hopf bifurcation encountered 
when 7 is increased from is the oscillation threshold: it 
is the minimal blowing pressure jth needed to initiate self- 
sustained oscillations (i.e. to play a note) when starting from 
zero and increasing quasi-statically the blowing pressure. 

3.4 Branch and continuation of Hopf bifur- 
cations 

This definition of a Hopf bifurcation can be written as an 
extended algebraic system G = where G reads: 

F(u, A,/i) 

Fu(u, A, -jcj(/) (10) 

A and /i are two parameters of interest (for instance the 
blowing pressure 7 and the reed opening parameter C), and 
(f) is the normed eigenvector associated with the purely imag- 
inary eigenvalue juj. 

Assuming that a given solution Xq = (-uq, ^o, ^0, Aq, /io) 
is a regular solution of G = 0, there exists a continuum of 
solutions X{jii) = (ii(/i), 0(/i), a;(/i), A(/i), /i) near Xq. Thus, 
the function X(/i) is a branch of Hopf bifurcations of our 
initial dynamical system. 

Using the same continuation techniques as in |3.2[ this 
branch of Hopf bifurcations can be computed. The reader 
is kindly referred to classical literature on the subject 
for proper theoretical definitions and other details about 
the continuation of Hopf bifurcations (see for instance 
Doedeli2Qj). 
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Figure 1: Impedance spectrum (reduced modulus) of a purely 
cylindrical bore of length L=57cm taking into account thermo- 
viscous losses and radiation and truncated to 18 modes as used 
in this study ( — ), analytical impedance spectrum without modal 
decomposition (-•-), impedance spectrum used in Wilson and 
Beavers [9 . ( ). 

4 Results 

In the current section, the method described in the previ- 
ous section is applied to the physical model presented in 
section [2] to investigate the oscillation threshold of a purely 
cylindrical clarinet. The resonator modal decomposition has 
been computed with the MOREESC[24j software, using 18 
modes. The corresponding input impedance spectrum is 
shown on figure [T] For comparison, the impedance spec- 
trum given by the analytical formulae of Backus [Ij and used 
in Wilson and Beavers [9j is also plotted. Important dif- 
ferences appear concerning frequency dependent damping 
(peaks heights) and harmonicity (peaks positions). 

4.1 Reed-bore interaction 

4.1.1 Comparison with previous results 

Previous works of Wilson and Beavers [9j and Silva et al.|12) 
showed how the resonance of the reed competes with the 
acoustical modes of the resonator for the existence of self- 
sustained oscillations: oscillation thresholds and emergent 
frequencies were then measured and numerically computed, 
for different values of the dimensionless product krL {kr = 
ujr/c is the corresponding wavenumber of the reed modal 
frequency, and L is the length of the bore -which is approx- 
imately the quarter of the first acoustical mode wavelength) 
called "tube parameter" in W&B. Notice that in these stud- 
ies, the parameter kr was kept constant and L was varied. 
These results showed that several regimes can be selected, 
depending on the product /c^L, only if the reed damping is 
very small. 

Figure [2] shows similar computations using our method, 
with the same model parameters as was used for the figure 




kL 

r 



Figure 2: Critical blowing pressure and frequencie of each Hopf 
bifurcation (labelled from 1 to 5) as a function of the "tube pa- 
rameter" krL (as used in W&B), or the frequency ratio Ur/c^i 
(secondary x-axis on top). Lower plot: critical blowing pressure 
7. Upper plot: frequency ratio = uj/ujr of the oscillations; 
passive resonances of the bore On = uJnjoJr are plotted for com- 
parison (-•-). Physical constants and parameters of the model 
were chosen according to fig. 1 in Silva et al.[12 : p=1.185g/L, 
c=346m/s, L=57cm, r=7mm, ^^=0.4, (^=0.13, Sr—^-, varying cjr- 
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Figure 3: Oscillation threshold : blowing pressure, regime selec- 
tion and frequency with respect to krL. The oscillation threshold, 
for a given abscissa, is given by the lowest blowing pressure of the 
five Hopf bifurcation branches (thick line). Lower plot: blowing 
pressure threshold 'jth- Upper plot: frequency ratio 0th = uJth/oJr 
of the oscillations. Complete branches of the five bifurcations are 
reminded (dashed lines). Same model parameters were used as 
fig. [2] A secondary x-axis on top of each plot gives the value of 
the frequency ratio uOr /ooi. 
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1 of Silva et al.[12j, only using a different impedance (as 
explained above). A secondary x-axis, on top of each plot, 
shows the corresponding ratio of the reed natural frequency 
uOr to the first resonance frequency of the bore uoi = Im(5i). 
Note that instead of varying L for a given cj^, the figure 
was computed by varying uj^ for a given L. The method 
advantage, here, is that no analytical development of the 
impedance spectrum is needed, unlike what was done by 
W&B to derive the characteristic equation. Thus, any other 
impedance could be used, for instance a measured one. 

Figure [3] illustrates how the oscillation threshold is de- 
duced from the previous figure: for a given abscissa, it cor- 
responds to the lowest of the five Hopf bifurcation branches, 
each corresponding to a register. In contrast with W&B, 
this figure clearly shows that even in the case of a strongly 
damped reed {q^ = 0.4), a regime selection occurs with vary- 
ing Wr. The term "strongly damped" was used by W&B for 
that value of but it seems to be a realistic value for a 
clarinet reed coupled to the player's lip (see van Walstijn 
and Avanzini[25j and Gazengel et al.[26j for numerical and 
experimental studies of the reed-lip system). A minimum of 
the threshold blowing pressure (lower plot) appears a little 
above ujr = cji, i.e. when the first air column resonance 
frequency is close to the reed one. The frequency of a given 
mode (upper plot) is close to the corresponding passive res- 
onance frequency ujn=lm{sn)^ for medium and high values 
of krL, but tends to the reed modal frequency for low values 
of krL. 



4.1.2 Relevance of the tube parameter krL 

While results of the previous figure are in good agreement 
with Silva's results (when thermoviscous losses are taken 
into account, as in figure 3 of Silva et al.[12j), several differ- 
ences lead to question the relevance of the representation: if 
krL is a characteristic parameter of the model, varying L for 
a given Ur should be equivalent to varying Ur with a fixed 
value of L. 

To answer this question, we recomputed the first regime 
for the same values of krL but with L=14.69cm (one then 
needs to recompute the modal coefficients of the input 
impedance). The results are plotted in figure |4j considering 
the upper plot, the frequency seems to be independent of L, 
as long as the product krL is kept constant (the plain line 
and dashed line are exactly superimposed); however, consid- 
ering the lower plot, the blowing pressure does not behave 
in the same way, for identical values of krL, depending on 
the value of L. Therefore, the tube parameter krL is not a 
characteristic invariant parameter of the model. 

One good reason for that is that in our impedance spec- 
trum, the peaks are inharmonic and have a frequency depen- 
dent quality factor and magnitude. Then varying L does not 
preserve the impedance peaks height and width, nor their 
spacing, which leads to a different balance in the competi- 
tion with the reed resonance. 
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Figure 4: Comparison of the first regime for L=57cm ( — ) and 

L=14.69cm ( ). Lower plot: dimensionless blowing pressure 

7 of the bifurcation. Upper plot: corresponding dimensionless 
frequency ratio = uj/ujr] the first passive resonance of the bore 
is reminded (•••)• 

4.2 Simultaneous influence of reed damp- 
ing and modal frequency 

In previous works [9 , 13j, figures similar to figure [2] were plot- 
ted for two cases, small and large values of the reed damping 
qr, which revealed very different behaviours. 

Because the method proposed in this paper leads to very 
short computation time, it is quite easy to loop the compu- 
tation for a series of qr values. For sufficiently close values 
of qr, this allows to draw a tridimensional plot in which 
the surface represents the critical value of the dimensionless 
blowing pressure 7 corresponding to a given Hopf bifurcation 
as a function of two parameters, the reed damping qr and 
the frequency ratio cj^/^i, as illustrated figurejsjfor the Hopf 
bifurcation of the first register. The plotted surface corre- 
sponds to a Hopf bifurcation locus when two parameters are 
varied (conversely to figure [2] where only one parameter is 
varied). The two plain, thick curves correspond to the limits 
of the domain considered: qr = 0.05 and qr = I. 

This figure allows to investigate the transition between 
large and light damping. It appears that for large values of 
ujr/oJi^ the critical blowing pressure is nearly independent 
of qr'. a closer look reveals a slight increase with qr. This 
is not surprising since Ur/uJi » 1 corresponds to the case 
where the first resonance frequency of the resonator is very 
small compared to uJr- Hence, only very large values of qr 
{qr » 1) would contribute significantly to an increase of 
the pressure threshold. However, for smaller values of ujr/oji, 
the value of qr becomes determinant. In the 2D manifold, 
there is a valley which becomes deeper when qr decreases. 
For a given g^, the bottom point of the valley (called 70 in 
Silva et al.[T2]) corresponds to the minimum threshold. It 
is plotted as a dot-dashed line on the plane (ujr/oji = 0), 



6 



0.5 




Figure 5: Surface giving the critical blowing pressure 7 of the 
Hopf bifurcation corresponding to the first register, with respect 
to the reed damping parameter qr and the reed modal frequency 
ujr (divided by the first bore resonance cji). The plain curves 
correspond to the limits of the chosen range for qr , i.e. qr = 0.05 
and qr = 1. The dot-dashed curve in the plane uJr/uJi = 
corresponds to the minimum 70, and the dot-dashed curve in 
the plane 7 = corresponds to the value of ujr/uoi at which it 



whereas the corresponding abscissa 



is plotted as a dot- 



dashed line in the plane 7 = 0. It shows that the minimum 
value 7o is an increasing function of qr . The same conclusion 

holds for ^ . 

Notice that the chosen range for the values of qr seems 
quite realistic, according to literature on the subject: in 
a numerical model, van Walstijn and Avanzini[25j reported 
parameters values equivalent to qr = 0.24 for one playing 
condition, whereas Gazengel et al.[26j found experimental 
values going from qr = 0.05 for the bare reed to qr = 1.54 
for a high lip pressure on the reed. Thus, the lower value 
qr = 0.05 corresponds to the limit case of a very resonant 
reed with no lip pressing on it, and the higher value qr = 1.00 
is high enough to cover a fairly good range of lip pressures. 

Now we illustrate that the results obtained on this model 
allow to estimate the range of validity of analytical formula 
obtained in approximated cases. For instance, in figure |6j 
the minimum 70 obtained through numerical continuation 
(without approximation), is compared with two analytical 
formula from Silva[13j corresponding to different approxi- 
mations: no losses in the cylindrical bore (here plotted with 
dotted line, corresponds to eq. (14) in Silva et al.[12j) and a 
single- mode resonator with losses (dot-dashed line, eq. (19) 
in Silva et al.|12j). 

It appears that only one mode with viscothermal losses 
leads to a more precise result than the undamped formu- 
lation with all the modes. However, the relative deviation 
from our results still reaches 10%. 



Figure 6: Minimum threshold 70 as a function of qr. Our nu- 
merical results ( — ), without approximation. Analytical results 
using two approximations: no loss single-mode resonator 

with losses (-•-). 



However, the oscillation threshold is not always given by 
the Hopf bifurcation corresponding to the first register. Con- 
sidering not only the first register, but the first five registers, 
leads to figure [7[a. Note that the number of registers is not 
arbitrary: at the right end of the figure, where ujr/oji 10 
{krL 16), the sixth resonance {ujq ^ llcji) is well beyond 
the reed resonance and thus cannot set up self-sustained os- 
cillations because its Hopf bifurcation is beyond 7 > 1. More 
registers should be considered if greater ujr values were to be 



explored. For sake of clarity, 70 and 



are not plotted. 



The plain lines in the plane {qr = 0.05) shows similar be- 
haviour as in figure [2] and has already been discussed. The 
other black plain lines depict a very different situation where 
the first register always correspond to the lowest 7. 

Figure [TJ-b shows the oscillation threshold, as defined by 
the first Hopf bifurcation encountered when increasing the 
blowing pressure 7 from 0, as a function of the reed damp- 
ing qr and the frequency ratio Ur/uJi. It is deduced from 
the previous one the same way figure [3] was deduced from 
figure [2] Intersections between the five surfaces result in an 
oscillation threshold surface with several local minima. 

It clearly shows how the reed damping, which can be con- 
trol by the player with his lower lip, plays a key role in the 
register selection, as previously reported by W&B[9j and 
Silva et al.[12j. It is also clearly visible that the range of qr 
for which a given register exists decreases with the index of 
this register. For instance, considering the fifth register, it 
can only be selected for qr < 0.3. 

The frequencies corresponding to the different registers 
are also calculated and plotted in figure |8] The frequency 
at threshold appears to be an increasing function of g^, but 
the influence of qr does not look significant in this 3D rep- 
resentation. When qr goes from to 1, the typical relative 



7 





Figure 7: a) Similar plot as in figure [s] for the first five registers: 
the critical blowing pressure 7 of each Hopf bifurcation is plotted 
with respect to the reed damping parameter qr and frequency 
ratio uJr/uJi- The plain curves correspond to the smallest and 
highest qr where a bifurcation occurs on the chosen range. 




Figure 7: b) The oscillation threshold extracted from figure |7]-a: 
blowing pressure at threshold 7^^, defined by the lowest of the 
five surfaces plotted in figure 7- a. 



Figure 8: Dimensionless frequency — uj/ujr at the bifurcation 
for the first five registers, with respect to the reed damping qr 
and the frequency ratio uJr /ooi. 

frequency deviations from the bore resonances are less than 
1%. However, a frequency shift more than 4% can be ob- 
served for the first register around uOr/ooi = 1.1. Such a ratio 
is quite unusual for a clarinet, but it is of great interest for 
reed organ-pipe manufacturer, where the reed natural fre- 
quency is close to the bore resonance, and the damping very 
small. 

4.3 Influence of the control parameter ( 

Let us remind here the definition of this parameter: ( = 
WZcy^2ho/ pKa- Its variations are mainly related to the 
changes of the reed opening parameter ho during the play. 
Those are directly driven by changes of the player's lip pres- 
sure and position on the reed. Thus it is a very important 
control parameter of the model. Notice that players modify 
both Qr and ( when changing the embouchure. 

Figure [9] shows the influence of ( on the blowing pressure 
and frequency at the oscillation threshold (expressed as the 
relative deviation to the corresponding bore resonance fre- 
quency). This control parameter happens to be critical for 
the regime selection : from very low values up to C=0.17, 
the flrst regime (with a frequency close to the flrst acousti- 
cal mode of the bore /o) is selected, while the fourth regime 
(/ 7/0) is selected for higher values. 

The control parameter ( also has noticeable influence on 
the frequency of the oscillations at threshold: whereas the 
first regime frequency deviation is less than 0.3% (5 cents) 
in the range of ( where it is the selected regime, the fourth 
regime frequency deviation is as high as 2.7% (46 cents) 
for C=0-8- This is a very important feature that has been 
highlighted in a paper by Guillemain et al.pT], where the 
lip stress on the mouthpiece of a saxophone were measured 
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Figure 9: Variations of the oscillation threshold with respect to 
(. Lower plot: critical, dimensionless blowing pressure 7; the low- 
est curve defines the oscillation threshold jth at a given abscissa. 
Upper plot: relative frequency deviation Auj/uJn = {uj — uJn) /oJn 
between the frequency at threshold ujth and the corresponding 
acoustical mode frequency ujn — Im(sn)- Model parameters: 
qr=OA, /r=1500Hz, S'^^Ocm^ 



while a player was playing, showing an adjustment of the 
lip stress on the reed in order to correct the tone shortly 
after the beginning of the oscillations. It could also explain 
the difficult reproducibility of measurements when fitting a 
clarinet or saxophone mouthpiece in an artificial mouth. 

The frequency deviations are monotoneous, decreasing 
functions of almost linear on the range of interest. Also, 
when C tends towards 0, the frequency at threshold tends 
to the corresponding bore passive resonance frequency. This 
result was to expect, since the boundary condition at the in- 
put end tends to a Neumann condition (infinite impedance). 
However, for very low values of the blowing pressure 
threshold ^th quickly increases, and eventually reaches 1 
which is the also the static closing threshold. In the case 
where '^th > 1, the reed channel is always closed and no 
sound is possible. 

4.4 Concurrent influence of Qr and ( 

Figure [To] shows the blowing pressure threshold with respect 
to the reed damping Qr and to the control parameter (. For 
a given pair (g^, ("), the oscillation threshold is given by the 
lowest surface among the five 2D manifolds corresponding 
to the different registers. 

Because the player's lower lip pressure and position on 
the reed, induce variations of both ( and Qr at the same 
time (the lip pressure being positively correlated with qr and 
negatively with it is interesting to follow the oscillation 
threshold along the horizontal plane ( = 1 — qr. On the 
one hand, for small values of and qr close to unity (i.e. 




Figure 10: Oscillation threshold with respect to the reed damp- 
ing parameter qr and the control parameter ^. The blowing pres- 
sure threshold 7th, is calculated the same way as in figure|7|b. The 
surface shading indicate the selected regime. 



a high lip pressure) the first regime is clearly selected; on 
the other hand, for a high value of ^ and small damping qr 
(i.e. a relaxed lip on the reed), it is the fifth regime that 
is clearly selected. In between, all regimes are successively 
selected except the second regime that is noticeably absent 
of the figure. 



4.5 Influence of the reed motion induced 
flow 

In the numerical results presented so far, the reed induced 
flow was not taken into account (5*^ = in all previous 
figures). However, as pointed out by Silva et al. [12j, it is 
important to take into account this additional flow, in order 
to predict accurately the emergent frequency at oscillation 
threshold. 



Figure [TT] shows the influence of the flow induced by the 
motion of the reed, through the value of the equivalent 
area Sr of the reed that participate to this flow. Previ- 
ous works by Dalmontp^, Kergomard[TT], or SilvafT ^ [TB ] 
report values around AI/=10mm, which corresponds, for a 
reed stiffness Ka=^.l^^Pa/m and a bore of radius r=7mm, 
to6'^=0.86cm2. Thus we scaled the variations of Sr between 
and 2cm^, using parameters values corresponding to a re- 
laxed embouchure, in order to well illustrate its influence. 

An important frequency deviation is shown as Sr in- 
creases, as well as a regime selection "cascade": the selected 
register is successively the fifth, the fourth, the third, and 
finally for high values of Sr the first. Finally, comparing 



with figure 10 (in the plane C=-75), the reed induced flow 
acts in a similar manner as the reed damping on the blowing 
pressure threshold and the regime selection. 
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Figure 11: Variations of the oscillation threshold with respect to 
Sr- critical blowing pressure 7 and relative frequency deviation 
from the bore resonances Auj/ujn for each of the five Hopf bifurca- 
tions. Parameters of the model: qr — 0.3, C = 0.75, fr — ISOOi/^. 

5 Conclusion 

In this paper, a clarinet physical model is investigated using 
numerical continuation tools. To the authors' knowledge, it 
is the first time that such method is used to compute the os- 
cillation threshold variations with respect to several model 
parameters. The reed dynamics and its induced flow are 
shown to have critical influence on the regime selection and 
the minimal blowing pressure necessary to bifurcate from 
the static regime and establish a steady-state periodic oscil- 
lation: a note. Previous works already gave useful insights 
concerning the influence of some model parameters on ease 
of play and intonation. The present work confirm and ex- 
tend these results to the case of a more complex model. 
Moreover, the method used allows to investigate variations 
of two parameters at the same time, instead of one, like in 
previous studies. 

In this approach, the parameters of the model are assumed 
to be constant or to undergo quasi-static variations. How- 
ever, in real situation, the player can modify some control 
parameters (e.g. blowing pressure, lip stress on the reed) 
at a timescale that might sometimes be comparable to the 
oscillations period. Guillemain[27j reported measured varia- 
tions of 7 on a time scale of a few milliseconds only, which is 
comparable to the period of an oscillation at Ib^Hz. How- 
ever, despite such a limitation, the results provided through 
numerical continuation are out of reach for direct time sim- 
ulations. 

Whereas the main results presented here concern a clar- 
inet model, it should be noted that the method itself is 
very general. No hypotheses (other than linear behaviour) 
is made on the resonator. Thus, other resonators can be 
studied by fitting the modal decomposition to its input 
impedance spectrum. For instance, extending to the case 



of the saxophone or taking into account the tone holes only 
requires to compute the corresponding modal decomposi- 
tion of the input impedance. Even the physical model can 
be modified. It only has to be written as a set of first or- 
der ordinary differential equations (and additional algebraic 
equations, if necessary). The method could also be used to 
compare models with each other. 

In comparison with the method of the characteristic equa- 
tion, the computations are much faster: Silva reported ten 
minutes of computation per branch, whereas it lasts only a 
few seconds in the present case. Moreover, the continuation 
algorithm used is very robust: strong variations (as varia- 
tions of 7 for low krL in fig. |2| do not require special care 
and are computed in a straightforward way. 

The same continuation method applied to the continua- 
tion of periodic solutions, as described by the authors in 
a conference paper [15j, also allows to compute the entire 
dynamic range of a given model of wind instrument, with- 
out any additional simplification, unlike previous works (see 
Dalmont[28j et al). From that perspective, the numerical 
continuation approach seems promising for the global inves- 
tigation of the behaviour of a given physical model of musical 
instrument. 

Applications of this work to instrument making are pos- 
sible: for instance, modifications of the geometry of the res- 
onator can be studied in terms of their influence on ease of 
play and intonation. Other applications concerning mapping 
strategies for sound synthesis are also of interest. Indeed, the 
estimation of the parameters of a model is a difficult task, 
especially when the parameters values should vary to repro- 
duce through sound synthesis typical behaviours of the in- 
strument modelled. On the one hand, direct measurements 
on a real player are most of the time highly complicated 
(for example, cj^, Qr or even when the lips of the player is 
pressed on the reed). On the other hand, inverse problems 
still appear to be limited to the estimation of parameter 
values from a synthesised sound. 

Thus, the approach presented in this paper offers the pos- 
sibility to know in advance the influence the parameter val- 
ues on some key features of the model behaviour around 
the oscillation threshold : ease of play, regime selection and 
fine intonation. Therefore, mapping strategies could be de- 
veloped for sound synthesis applications. They would con- 
sist in binding different parameters in order to let a player 
modify one of them while maintaining the same playing fre- 
quency (at threshold) for a given note for example, or while 
preserving the ease of play on a whole register. 
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